Dynamic instabilities of fracture under biaxial strain using a phase field model 
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We present a phase field model of the propagation of fracture under plane strain. This model, 
based on simple physical considerations, is able to accurately reproduce the different behavior of 
cracks (the principle of local symmetry, the Griffith and Irwin criteria, and mode-I branching). In 
addition, we test our model against recent experimental findings showing the presence of oscillating 
cracks under bi-axial load. Our model again reproduces well observed supercritical Hopf bifurcation, 
and is therefore the first simulation which does so. 



5h 
43 



. 

C ■ 

o ■ 
> 

m ■ 

: 

in . 
(N ■ 

o ■ 

: 

o . 
-i— > ■ 

a ■ 

s; 

T3 ■ 

c : 
o . 

o ■ 



X 



In recent years, the physics community has seen a 
rebirth of interest in the problem of dynamic fracture. 
This rebirth was kindled by a series of experiments re- 
vealing that the current engineering approach to crack 
propagation, namely the coupling of linear elasticity 
to an empirical energy balance law for crack tip mo- 
tion, cannot account for the richness of actual fracture 
phenomenology^. Specifically, dynamical instabilities 
which drive the system away from a single crack propa- 
gating in a straight line require more sophisticated atten- 
tion to the actual tip region, the so-called process zone. 

Given the above, it is clear that one needs a frame- 
work which can couple local degrees of freedom involved 
in breaking inter-atomic bonds to global elasticity. One 
approach, referred to as phase-field modeling of fracture, 
accomplishes this task by introducing an order-parameter 
field (the degree of "broken-ness" ) which then couples 
to the elastic strain in a manifestly continuum-level for- 
mulation. The fact that the system does not need to 
be placed on a lattice avoids dynamical artifacts asso- 
ciated with the breaking of translational and rotational 
symmetry 0. One such phase-field model due to Karma, 
Kessler and Levine (KKL) 3] has been shown to correctly 
encompass much of the expected behavior of mode III 
(out-of-plane) cracks Q. 

Here, we extend the KKL model to full vector elas- 
ticity and test its genericity. Our major interest is in 
seeing whether a recently discovered supercritical Hopf 
bifurcation to oscillating cracks under biaxial loading [5j 
(see Fig. is in fact reproduced by KKL; we will see 
that in fact it is. This first-ever successful simulation of 
the crack oscillation has the dual benefit of demonstrat- 
ing that the instability is not dependent on any special 
properties of the specific materials used in the experi- 
ment (rubber) and also of giving us more confidence in 
the phase-field methodology. 

We start with the KKL phase field model 3]. Here, a 
sheet of fractured elastic material is represented by the 



and u y derives from a modified elastic energy: 



elastic displacement field u x 



and by a phase field 



variable <p that can be interpreted as the proportion of 
intact inter-atomic links. The evolution equation of u x 
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where g = (4 — 3<fi)4> 3 is a function of (j> chosen such that 
3(0) = 0, 0(1) = 1, and g'(0) = g'(l) = 0; this specific 
choice is discussed in The tensor is the strain 
tensor which has the following form: 
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The evolution equations for u x and u v are then: 
pduUi = - — 

OUi 

The corresponding evolution equation of <f> is 



r^ = A^-^M-«(^-e c ) (4) 



where V(<p) — 4(0(1 — </>)) 2 is a double well potential and 
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where KL ame = ( A + /j)/2 is the modulus of compression 
for a plane strain configuration; a is an arbitrary coef- 
ficient chosen to be bigger than 1 . This breaking of the 
symmetry between compression and extension is a key in- 
gredient not present in previous phase field approaches to 
in-plane fracture@,0- In our model, If a material is sim- 
ply compressed, having a > 1 will guarantee that it will 
not break; also, compression will increase the threshold 
needed for inducing a fracture through shear. The pa- 
rameters values used here are: e c = l,r = 5,A = /i = l 
and a = 1.5. 

We proceed to study this model computationally in a 
box of width W and length L. Simulations were per- 
formed using a grid spacing of Sx = 0.15 and a time step 
of St = 0.001 ( Decreasing to Sx = 0.075 and St = 0.0005 
leads to no significant difference in the results). The time 
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FIG. 1: Left: Definition of in-plane loading modes: top - 
mode II loading, bottom - mode I loading. Right: Geometry 
of the experiment reported n |f| . An elastic sheet is extended 
in both x and y directions, such that the local strain is bigger 
in the y direction than in the x direction. The crack speed 
is denoted by c. In simulations performed here W was taken 
equal to 60 and boundary conditions at y = and y — W 
were chosen so that u x was kept constant in time and hence 
so that a vertically propagating crack can not reach y = W 
or y = 0. 

stepping scheme was the forward Euler method while the 
spatial operators were computed using a discretization 
that conserves the discretized energy of Eq. ^ (The use 
of other scheme lead to long term numerical instabilities). 
In our simulations, we used both fixed grids of different 
sizes along the x axis and a grid moving with the fracture 
tip along the x axis. Boundary conditions for the fixed 
grid were as follows: at the y — W (y = 0) boundary, 
u y was kept equal to A y (0) and on both lateral edges 
u x = xA x /L. At the x = L (x = 0) boundary, u x was 
kept equal to (0) and no flux boundaries were used for 
the u v field. In the case of the moving grid, the boundary 
conditions at y = W and y — were unchanged whereas 
the boundary conditions at the horizontal ends of the grid 
were modified. First we introduced an artificial viscosity 
in a thin layer at both ends (This mimics partially ab- 
sorbing boundary conditions). Also, at the leftmost end 
the u x field was kept constant between each displacement 
of the grid. We checked that those modifications did not 
affect the behavior of the crack when compared to results 
obtained using a long enough fixed grid, thereby allow- 
ing us to simulate the crack propagation along an infinite 
strip along the x axis. The initial conditions were con- 
structed as follows: a small initial crack was created by 
setting <f> — in a small region of fixed width and vari- 
able length and by letting the system evolve following a 
damped version (see later) of the evolution equation of 
elasticity (while <p was kept constant) until a stationary 
state was reached. 

We first tested the model for damped dynamics ob- 
tained by replacing <9 tt by dt- For pure mode I load- 
ing, the fracture began to propagate once the imposed 
elastic energy was higher than the fracture energy (see 
later). In addition, in this case, the stress intensity fac- 
tor at the tip of a steady crack was found to be constant 
for various loading configuration, hence obeying the Ir- 
win criterion|9j (data not shown). More interestingly, 




FIG. 2: Left: contour plot of the line (f> — 0.5 taken at 1. 
time unit (tu) intervals during the branching of a crack under 
mode I loading with A y = 18. and J3 = 1. The angle of the 
two branches at the branching site is approximately 40° and 
the critical speed for which the instability occurs is approxi- 
mately 0.5 while the shear speed is 1, which is in agreement 
with predictions in During the evolution of the system, 
the lower branch will recede while the upper branch will prop- 
agate and branch irregularly. The simulations are performed 
using the undamped model. Right: contour plots of (f> = 0.5 
taken at 40. time units (t.u.) during the propagation of a 
damped crack in a medium where pure mode II loading is 
applied. The initial crack is straight and oriented along the 
x axis. The crack propagates in a direction that nullifies the 
mode II stress intensity factor. Inset: same parameters with 
a — 0. One can see, that in this case, the crack, branches and 
propagates in two directions symmetrical with respect to the 
x axis. The second branch propagates in a region where the 
material is strongly compressed. Simulations are performed 
using the overdamped model. 

numerical simulations in the case of pure mode II load- 
ing showed that the model respects the local symmetry 
principle: the fracture propagates in the direction which 
nullifies the stress intensity factor for mode II (see Fig. 
2b). Hence this model reproduces well the behavior of a 
single crack in the damped regime. Note that this result 
depends on our asymmetry parameter a; allowing break- 
age under compression leads to a model which does not 
follow the local symmetry principle (Fig. 2b, inset). The 
results for different a ranging from 1 to 2 were similar. 

Next, we turn briefly to results obtained in the dy- 
namic (non-damped) case under pure mode I loading. As 
expected by analogy with the results of 3] for the case 
of mode III loading, the initiation of crack propagation 
appeared at the Griffith threshold with good accuracy. 
Indeed, according to the Griffith criterion, one would ex- 
pect a crack to begin to propagate for a value of the y 
extension A y bigger than A c which is solution of: 

W ((^ + v)^j=2^ dcj> ^2(1 - g(4>)) + V(<j>) 

(6) 

This formula was derived in |3( and follows from the 
asymptotic solution of the model in the region far behind 
the crack tip. For the parameter values used here, this 
threshold is at 9.5 whereas in our simulations the crack 
begins to propagate for A y bigger than 9.7±0.1; this 2% 
discrepancy is due to discretization and finite width ef- 
fects and gives some measure of the accuracy of our com- 
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FIG. 3: Fracture tip trajectories for A y — 10 and r = 
WA X /L — 8.0 during the transient regime (a) and when a sta- 
tionary state is reached (b). The steady-state tip trajectory 
matches a simple sinusoid, (c) Tip trajectory for A y — 10. 
and r = WA X /L = 7.25, exhibiting/ damped oscillations. In 
(b) the mean tip speed is 0.362 and the mean tip speed along 
the x axis is 0.354. In (c), the tip speed is 0.347 once straight 
propagation has resumed. The transverse and longitudinal 
speeds of sound are respectively 1. and 1.22. The period of 
the oscillations in (a) is about 831 t.u and the corresponding 
wavelength is 294 s.u 

putations. The speed of the stable crack behaved quali- 
tatively as expected when loading was increased. When 
loading was further increased, the dynamic fracture ex- 
hibits a branching instability and a secondary crack be- 
gins to propagate (see FigEJ). This branching instability 
is compatible with what has been observed experimen- 
tally in a wide range of materials [lfj, 111) . A complete 
analysis of the branching phenomenon will be addressed 
in further work. 

We now turn to our major interest here, the case of 
a dynamic crack propagating under biaxial stress. Ex- 
perimental work by Deegan et al. [j| has shown that for a 
given imposed strain in the y direction (see Fig. there 
is a threshold value of the x strain for which the crack 
propagation is no longer straight; instead, the crack tip 
position begins to oscillate. In fact, the instability ap- 
pears to be supercritical and the tip trajectory is well- 
approximated by a sinusoidal line with finite wavelength 
and amplitude. Recall that in our calculations strains 
are applied by moving the rightmost border of the sheet 
by A x and by moving the top border by A^; hence, if 
A x = 0, the system is set to pure mode I. The experi- 




FIG. 4: Elastic energy landscape during the propagation of 
an oscillating crack. The white region corresponds to the 
broken region (<j> < 0.5) and the gray scale represents the 
density of elastic energy with dark regions corresponding to 
high elastic energy density. Parameter values are as in figure 
0(a). Between the arches of the sinusoidal line of the crack, 
the elastic energy density is lower than in the vicinity of the 
upper and lower boundaries of the slice. 

mental results translate into the prediction of a Hopf bi- 
furcation that should occur as we cross a threshold value 
of r = WA X /L, with A y being fixed. 

The results of our numerical simulations for two sets 
of A, /i and different values of A y faithfully reproduces 
the aforementioned phenomenology. Namely, the frac- 
ture tip trajectory indeed undergoes a Hopf a bifurcation 
when the x extension is increased over a threshold value 
that depends on both parameter regime and the vertical 
extension. This bifurcation is characterized by the fact 
that below threshold, the tip position shows damped os- 
cillations (see figure [3J| and ends up propagating along a 
straight line, whereas above threshold those oscillations 
are amplified and the restabilized state corresponds to 
the situation where the fracture tip oscillates at a finite 
wavelength with a finite amplitude (see figure 01 . We 
checked that this instability was not due to waves reflect- 
ing at the boundaries that can create periodic markings 
called Wallner lines^^]. Indeed, the expected wavelength 
of such markings would be A = Wv/cy/1 — v 2 /c 2 , that 
is about 10 s.u. while the wavelength observed here is 
about 300 s.u. This is confirmed by the fact that switch- 
ing to quasi-absorbing boundary conditions at the y = W 
and y = lines does not affect the oscillations. 

When the restabilized state is reached, the trajectory 
of the tip is almost indistinguishable from a sinusoidal 
line, as in One can also note that the horizontal 
tip speed oscillates with a frequency equal to two times 
the frequency of the vertical position, so that the max- 
imum of the horizontal tip speed are reached when the 
instantaneous tip velocity is directed along the x axis. In 
addition, the tip speed tangent to its trajectory is kept 
almost constant (up to numerical errors) and for different 
values of r (r varying between 7 and 8), we did not find 
significant changes in the tip speed (less than a few %). 
A picture of this state is presented in Fig. 4. 

We now describe the changes in the oscillating restabi- 
lized state when the strain along the x axis is increased. 
As seen in figure the amplitude of the oscillations be- 
haves like \/ A x — A xc close to threshold, which is con- 
sistent with a supercritical Hopf bifurcation, as experi- 
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FIG. 5: (a): Square of the amplitude of the oscillations as 
a function of r = WA X /L. (b): frequency ui (+) of the os- 
cillations. The wavelength (not shown) varies like the inverse 
of uj since the crack speed is not varying much (about 10~ 3 ) 
over the parameter range presented here. 



mentally observed in |£g. The wavelength and period of 
the tip oscillation decreases slightly when — A xc is 
increased; this differs from results in [{| where the wave- 
length increases when A x — A xc . This may be due to 
non-linear elasticity effects present in the experiment. 

We also performed numerical simulations with biax- 
ial strain and a = O.Thc results obtained did not differ 
significantly from that observed for a = 1.5. This is 
somewhat surprising, since the simplest interpretation of 
the oscillations observed here suggests that the mecha- 
nism is at least partially similar to the one underlying 
oscillations of a quasi-static crack propagating in a ther- 
mal gradient |l3l Il4| . In that case, theoretical work has 
explained the transition to an oscillatingcrack using a 
(modified) principle of local symmetry |l5l llfij : chang- 
ing the tip direction induces a mode II component which 
causes further deviation from the original line. We have 
already shown that this principle does not apply for a 
symmetric model for pure mode-II loading. Perhaps the 
explicit breaking of the symmetry by the mode-I part of 
the driving is enough to suppress the unphysical com- 
pressional breaking for the case of a — 0, and hence this 
model still exhibits the Hopf bifurcation. In support of 
this, we verified that even in this case, a crack tip with 
damped dynamics will obey the principle of local sym- 
metry, if in addition to a pure mode II load one adds a 
small mode I extension. We should note, though, that we 
never observe oscillations with damped dynamics, even 
when r was set to very close to A y , i. e. close to hy- 
drostatic strain. However increasing r (up to 50), i.e. 
increasing the dissipation at the crack tip did not affect 
significantly either the instability wavelength of the os- 
cillating crack or the threshold but did reduce the crack 
speed (from 0.362 to 0.05 (r = 50)). Also, changes in A y 
(A y = 12, 14, with t — 20 to avoid branching) did not af- 
fect significantly the wavelength of the oscillating crack 
but did change the threshold. Interestingly, the wave- 
length scales linearly with W . Hence, it seems that the 
saturation of the amplitude is governed by the interaction 
of the tip with the sidewall. Underlying the actual insta- 



bility is perhaps the simple fact that an oscillating crack 
will alleviate extra elastic energy under biaxial strain. 

In summary, this paper shows that the extension of 
the KKL phase field model of crack propagation to full 
vector elasticity qualitatively reproduces the different in- 
stabilities observed when considering the propagation of 
cracks. One should note that with an extremely simple 
model based on generic physical considerations we were 
nonetheless able to reproduce the variety of observed pat- 
terns. This lends us confidence in the entire modeling 
approach and suggests that we proceed in two comple- 
mentary directions. First, we can continue to investi- 
gate the phenomenology of KKL, specifically looking at 
the interaction of different cracks (the phase field method 
can easily deal with intersecting interfaces) and also truly 
three-dimensional effects^! ■ At the same time, it is time 
to begin understanding how to combine this method with 
microscopic interaction data about specific materials so 
as to enable the building of more quantitatively reliable 
models of dynamic fracture. 
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